function HCP_KRR_predictable_behaviors(KRR_dir, maxKRR_iter, Nperm, test_metric, intrim_csv, outmat, ...
    restricted_csv, subj_ls, bhvr_ls, colloq_ls)

    % HCP_KRR_predictable_behaviors(KRR_dir, maxKRR_iter, Nperm, test_metric, intrim_csv, outmat, ...
    %    restricted_csv, subj_ls, bhvr_ls, colloq_ls)
    %
    % Permute behavioral scores across subjects, considering the family structures (multi-level block 
    % permutation). Muilti-level block permutation requires the FSL PALM package.
    % Run kernel ridge regression (KRR) on the permuted scores to test the predictability of the original 
    % KRR model.
    %
    % Inputs:
    %   - KRR_dir
    %     The directory which contains the KRR trained models and testing results.
    %   - maxKRR_iter
    %     Maximal random seed used to split the training-test folds for performing KRR, e.g. 400.
    %   - Nperm
    %     Number of premutations.
    %   - test_metric
    %     The accuracy metric used to perform statistical testing. Choose from 'predictive_COD' and 'corr'.
    %   - intrim_csv
    %     Full path of the intermediate csv file generated by hcp2block2 function.
    %   - outmat
    %     Full path of the output mat file storing the behaviors whose actual prediction accuracy was significantly
    %     above chance.
    %   - restricted_csv 
    %     Full path of the HCP restricted csv file, containing the family structure information. 
    %   - subj_ls
    %     Full subject list (absolute path).
    %   - bhvr_ls 
    %     List of behavioral measures for which matched AA and matched WA can be found (absolute path).
    %   - colloq_ls
    %     List of colloquial names of behavioral variables. The colloquial names should correspond to the
    %     behavioral names in "bhvr_ls".
    
    addpath(fullfile(getenv('CBIG_CODE_DIR'), 'external_packages', 'matlab', 'non_default_packages', 'palm', 'palm-alpha109'))
    
    iter = 40;
    
    %% load shared files
    [bhvr_nm, nbhvr] = CBIG_text2cell(bhvr_ls);
    [colloq_nm, ncolloq] = CBIG_text2cell(colloq_ls);
    if(ncolloq ~= nbhvr)
        error('Number of behavioral names is not equal to number of colloquial names.')
    end
    [subjects, nsub] = CBIG_text2cell(subj_ls);
    alpha_FDR = 0.05;
    
    %% create permutations
    % 1. use `hcp2blocks.m` to generate block definitions
    % 2. use `palm_tree.m` to create exchangeable tree
    % 3. use `palm_permtree.m` to generate permutation indices
    idx_perm_out = fullfile(KRR_dir, 'Pset.mat');
    if(~exist(idx_perm_out, 'file'))
        B = hcp2blocks(restricted_csv, intrim_csv, false, dlmread(subj_ls));
        Ptree = palm_tree(B);
        Pset = palm_permtree(Ptree, Nperm+1, false, true);
        save(idx_perm_out, 'Pset', 'B')
    else
        load(idx_perm_out)
    end
    
    %% load data for KRR; calculate permuted KRR accuracies
    metrics = {'corr','COD','predictive_COD','MAE','MAE_norm','MSE','MSE_norm'};
    parfor b = 1:nbhvr
        KRR_perm(b, KRR_dir, maxKRR_iter, Nperm, Pset, bhvr_nm, metrics)
    end
    
    %% calculate p value for each behavior
    p_perm = zeros(nbhvr, 1);
    avg_stats = [];
    avg_null_stats = [];
    for b = 1:nbhvr
        curr_avg_stats = [];
        curr_avg_null_stats = [];
        for i = 1:maxKRR_iter
            opt_fname = fullfile(KRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, ...
                ['final_result_' bhvr_nm{b} '.mat']);
            if(~exist(opt_fname, 'file'))
                continue
            end
            opt = load(opt_fname);
            acc_out = fullfile(KRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'perm.mat');
            load(acc_out)
    
            curr_avg_stats = cat(1, curr_avg_stats, mean(opt.optimal_stats.(test_metric), 1));
            curr_avg_null_stats = cat(1, curr_avg_null_stats, mean(stats_perm.(test_metric), 1));
        end
    
        avg_stats = cat(2, avg_stats, curr_avg_stats);
        avg_null_stats = cat(3, avg_null_stats, curr_avg_null_stats);
        p_perm(b) = length(find( mean(curr_avg_null_stats,1) > mean(curr_avg_stats,1) )) ./ Nperm;
    end
    
    %% multiple comparisons correction
    H = FDR(p_perm, alpha_FDR);
    sig_perm_idx = sort(H);
    sig_behaviors = bhvr_nm(sig_perm_idx);
    
    switch test_metric
    case 'predictive_COD'
        sig_COD = sig_perm_idx;
        p_COD_perm = p_perm;
        save(outmat, 'p_COD_perm', 'H', 'sig_COD', 'sig_behaviors', 'avg_stats', 'avg_null_stats');
    case 'corr'
        sig_corr = sig_perm_idx;
        p_corr_perm = p_perm;
        save(outmat, 'p_corr_perm', 'H', 'sig_corr', 'sig_behaviors', 'avg_stats', 'avg_null_stats');
    end
    
    
    rmpath(fullfile(getenv('CBIG_CODE_DIR'), 'external_packages', 'matlab', 'non_default_packages', 'palm', 'palm-alpha109'))
    
end
    
    
    
    
    
    
function KRR_perm(b, KRR_dir, maxKRR_iter, Nperm, Pset, bhvr_nm, metrics)
    
    for i = 1:maxKRR_iter
        opt_fname = fullfile(KRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, ...
            ['final_result_' bhvr_nm{b} '.mat']);
        if(~exist(opt_fname, 'file'))
            continue
        end
        opt = load(opt_fname);

        flag = exist(fullfile(KRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'FSM'), 'dir');
        if(flag)
            load(fullfile(KRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'FSM', 'FSM_corr.mat'));
        end

        load(fullfile(KRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, ...
            ['no_relative_10_fold_sub_list_' bhvr_nm{b} '.mat']));
        Nfolds = length(sub_fold);
        
        acc_out = fullfile(KRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'perm.mat');
        if(exist(acc_out, 'file'))
            continue
        end
        
        for f = 1:Nfolds
            y_reg = load(fullfile(KRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'y', ...
                ['fold_' num2str(f)], ['y_regress_' bhvr_nm{b} '.mat']));
            opt_lambda = opt.optimal_lambda(f);
            
            test_ind = sub_fold(f).fold_index==1;
            train_ind = ~test_ind;
            N_train = length(find(train_ind));
            N_test = length(find(test_ind));
            
            if(flag)
                K_train = FSM(train_ind, train_ind);
                K_test = FSM(test_ind, train_ind);
            else
                load(fullfile(KRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'FSM_innerloop', ...
                    ['fold_' num2str(f)], 'FSM_corr.mat'))
                K_train = FSM;
                clear FSM
                load(fullfile(KRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'FSM_test', ...
                    ['fold_' num2str(f)], 'FSM_corr.mat'))
                K_test = FSM(test_ind, train_ind);
                clear FSM
            end
            
            % compute the part of parameters that are not dependent on y
            % so that they are be shared across all permutations
            K_r = K_train + opt_lambda*eye(N_train);
            X = ones(N_train,1);
            inv_K_r = inv(K_r);
            beta_stable = (X' * (inv_K_r * X)) \ X' * inv_K_r;
            
            % for each permutation, calculate prediction accuracies
            for p = 2:Nperm+1
                y_perm = y_reg.y_resid(Pset(:,p));
                y_train = y_perm(train_ind);
                y_test = y_perm(test_ind);
                
                beta = beta_stable * y_train;
                alpha = inv_K_r * (y_train - X * beta);
                
                y_p = K_test * alpha + ones(N_test,1) .* beta;
                for k = 1:length(metrics)
                    stats_perm.(metrics{k})(f,p-1) = ...
                        CBIG_compute_prediction_acc_and_loss(y_p, y_test, metrics{k}, y_train);
                end
            end
        end
        save(acc_out, 'stats_perm')
    end
    
end